Stripe formation in bacterial systems with density-suppressed motility 
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Engineered bacteria in which motility is reduced by local cell density generate periodic stripes of 
high and low density when spotted on agar plates. We study theoretically the origin and mechanism 
of this process in a kinetic model that includes growth and density-suppressed motility of the cells. 
The spreading of a region of immotile cells into an initially cell-free region is analyzed. From the 
calculated front profile we provide an analytic ansatz to determine the phase boundary between 
the stripe and the no-stripe phases. The influence of various parameters on the phase boundary is 
discussed. 

PACS numbers: 87.18.Hf, 87.23.Cc 
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Biological systems exhibit a wide variety of exquisite 
spatial and temporal patterns. These patterns often 
play vital roles in embryogenesis and development [l|, [^j . 
In addition, colonies of bacteria and simple cukaryotes 
also generate complex shapes and patterns Typi- 
cally, these patterns are the outcome of coordinated cell 
growth, movement, and differentiation that involve the 
detection and processing of extracellular cues Q • 

These experimental observations have triggered exten- 
sive mathematical modeling. A large body of theoretical 
work is devoted to pattern formation by chemotactic bac- 
teria. On the mean-field level, these phenomena can be 
described by Keller-Segel type reaction-diffusion models 
In many instances, the models invoke non-linear 
diffusion of the cells where the diffusion coefficient in- 
creases with the local cell density 0, [l^l . 

Recently, it was theoretically proposed that the oppo- 
site case of density suppressing motility could also lead 
to patterns via a "self-trapping" mechanism [13, Q. In 
parallel, we have explored such a system experimentally, 
using a synthetic biology approach [l5|. The density- 



suppressed motility was introduced into the bacterium 
E. coll by having it excrete a small (and rapidly de- 
graded) signaling molecule AHL, such that at low AHL 
levels, these cells perform random walks via their swim- 
and-tumble motion 16[ and are "motile" , while at high 
AHL levels, these cells tumble incessantly, resulting in a 
vanishing macroscopic motility and becoming "immotile" 
[Fig. Ha)]. 

On agar plates, these engineered bacteria form highly 
regular and stable stripe patterns consisting of periodi- 
cally alternating regions of high and low cellular densities 
[Fig. [TJb)]. A thorough characterization of these spatial 
patterns gave rise to the following key experimental ob- 
servations (i) Regulation of cell motility by AHL 
is essential for pattern formation; (ii) Cells are motile at 
low densities and immotile at high densities; (iii) Bacteria 



form stripes sequentially in one- and two-dimensional ge- 
ometries when expanding into an initially cell-free region; 

(iv) Random initial conditions do not give rise to stripes; 

(v) Chemotaxis is not required for pattern formation; (vi) 
The stripe patterns depend on the magnitude of the un- 
repressed cellular motility in the low density limit: Upon 
decreasing this magnitude the system makes a transition 
from a phase with spatially periodic stripes (the stripe 
phase) to the no-stripe phase, through a region with a 
finite number of stripes. 

As demonstrated in all the experimental obser- 
vations can be reproduced by a three-component model 
that (i) describes the cellular motion as random walk 
with an abrupt AHL-dependent motility coefficient, (ii) 
takes into account the synthesis, diffusion, and turnover 
of AHL, and (iii) implements the consumption and diffu- 
sion of the nutrient due to cell growth and the limitation 
of growth in the absence of nutrient. 

Despite the success of this model, the origin and mech- 
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FIG. 1: (a) The engineered bacterium cells execute "random 
walks" at low densities but become immotile at high densi- 
ties, (b) This coupling between density and motility leads to 
the formation of stripes with periodic density variations on 
agar plates [31 . Initial cell seeding was done (at the position 
indicated by the arrow) 30hr before the picture was taken. 
Bar corresponds to a length of 5mm. 
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anism of the pattern formation process remain unclear. 
In this paper, we describe a simphfied two-component 
model to study the essential features of stripe formation 
analytically. In terms of the concentration h{x, t) of AHL 
and the cell density p{x, t) at position x and time t, the 
dynamical equations are given by, 
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The first equation describes production (with rate a), dif- 
fusion (with diffusion coefficient Dh) and turnover (with 
rate /3) of AHL. It is clear that in spatially homogeneous 
situations, /i cx p in the steady state, hence the name 
"quorum sensor" for AHL. The second term on the right 
hand side of Eq. © describes logistic bacterial growth at 
rate 7 and with a saturation density ps. The reduced 
growth rate at high densities approximates the nutri- 
ent depletion effect in the experiments. The stochastic 
swim-and-tumble motion of the bacteria is modeled as a 
diffusion-like term on the right hand side of Eq. ([2]) . The 
experimentally measured values of all parameters can be 
found in Ref. (l5| . 

The motility function nih) explicitly depends on h. It 
takes into account the repressive effect of AHL concen- 
tration (and hence cell density) on cell motility. The 
interaction term in Eq. ^ can be obtained by either 
generalizing the coarse-graining procedure of Ref. 13| or 
adopting the master equation approach of Ref. [l7j to an 
/i-dependent motility. In fact, such an analysis yields a 
mixture of two terms dx[p(h)dxp) and dx{pdxp{h)) (for 
details see SM) . For simplicity we focus on the above cou- 
pling, but our main conclusions are not affected by this 
(for details see SM). 

Measurements of bacterial diffusion at the population 
level show that p drops abruptly from a value Dp to 
Dpfi <^ Dp as h increases beyond a threshold ho. As 
simulation results of Ref. [l3| did not depend sensitively 
on the value Dp o, we shall set Dp o = 0. Thus, we con- 
sider the form p{h) = Dp for h < ho — w and p{h) = 
for h > ho with a linear decrease of p for the transition 
region ho ~ w < h < ho with hp ^ w ^ 0. 

As demonstrated in Ref. [15|, this two-component 
model is able to initiate stripe patterns in a growing 
bacteria colony and maintain them for a while; but the 
stripes are eventually lost after long times when cell den- 
sities reach ps throughout the system. The latter behav- 
ior deviates from the experimental system where stripes 
are frozen in upon nutrient exhaustion. Nevertheless, the 
model correctly captures the dynamics at the propagat- 
ing front where new stripes are formed. The simplic- 
ity gained enables analytic treatment that clarifies con- 
ditions for spontaneous stripe formation in the system. 

Consider a one-dimensional bacterial colony develop- 
ment as depicted in Fig. [TJb). Initially, cell density is 



low on the plate and all cells grow and freely diffuse. 
As growth proceeds cells at the center aggregate. The 
increased cell population boosts the local AHL concen- 
tration, driving it eventually above ho so that cells inside 
the aggregate become immotile. At the same time, this 
high density region expands outward by absorbing cells 
moving from surrounding low-density regions into the ag- 
gregate. Depending on the parameter values of the sys- 
tem, the high density region expands either stabl y a s a 
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front or exhibits instability that results in stripes 

We now take a closer look at the low-density region 
that precedes the advancing aggregate, whose cell den- 
sity profile is calculated later (see Fig. [5]). The size of 
this motile cell population is maintained by a dynamic 
balance between cell growth within and loss to the ag- 
gregate in the contact zone. Due to absorption by the 
aggregate, cell number is low in the contact region. By 
virtue of Eq. ([2]) , the maximum density p,„ of motile cells 
is found at a distance Lp = yjDpj^ from the aggregate, 
while the cell diffusion flux into the aggregate is given 
by J ~ DpPrajLp. Meanwhile, the expansion speed c of 
the aggregate satisfies J = cpc where pc is the density 
drop across the aggregate boundary. Hence quite gener- 
ally pm — Pci i-e., the cell density profile in the motile 
region scales with the density at the edge of the aggre- 
gate where the AHL concentration is at the threshold 
value ho- A quantitative calculation is then required to 
determine whether the AHL concentration rises to the 
threshold again at pm- As we shall see below, the answer 
depends on how the diffusion length = ^ D^j fi of 
AHL molecules (i.e., the typical distance travelled by an 
AHL molecule before degradation) compares to Lp. 

We will analyze a rescaled version of the model ([T])- 
([2]) that only depends on dimensionless quantities. We 
measure length in units of Lp, time in units of I/7, p in 
units of pho/oi and h in units of ho. All dimensionless 
quantities are denoted by a hat (e.g. t = etc.). For a 
steadily propagating front at speed c, the density profiles 
p and h become functions of z = x — ct. We set the front 
position at z = such that cells are immotile for z < 
(region I) and motile for z > (region II) . 

The cell density profile in region I is easily obtained 
by integrating Eq. ^ with the boundary condition 
/5i(-oo) = Ps = aps/iPho), 
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where pc = Pi{^ ) is the scaled cell density at the edge 
of the aggregate. The experimental system of Ref. [15[ 
has a ~ 4. 

In region II, the marginal stability criterion [isj yields 
P\\{z) ^ for z ':$> 1 with the selected wave speed 
c — 2. With this choice, Eq. ^ takes on the following 
form (except within a distance w from the interface). 



Pn + 2pii + pii(l-pii/p,) = 0, 



(4) 



where the prime denotes d/dz. For the form of 
described, one finds /5ii(0"'") = as w — J- (see SM). 
Matching the diffusional flux from the motile side with 
the speed of the immotile front yields the second condi- 
tion at the interface: /5{j(z = O"*") = 2pc- Thus, pii{z) is 
a non-nionotonic function, rising for small z before de- 
caying exponentially for large z. 

The AHL profile is determined from the cell density 
profile as 

/oo 
dzip{zi)Gh{z - zi), (5) 
-OO 

with the Green's function, 

Gh{z) = -(1 + A.A)e-^/^''e(^+^'*^'l^l/^'V2, (6) 

and A =-[! + (! + -D;,/3)i/2]/i)^. 

Due to the nonlinearity in Eq. (jl]), an exact solution for 
pii is not possible and we shall analyze the problem in an 
expansion in e = Pc/Ps = Pel Ps- We shall first consider 
the limit e — ?■ 0. The solution (|3]) is then approximated 
by pi(z) ~ = Pce"^/^. The linear form of Eq. (g]) 

together with the matching conditions at z = yields 
/5|™(z) = 2pcze~^. Inserting into ([5]), we obtain 

the AHL profile to the zeroth order in e, 

}^{z) = pp(^—^e-'- + he-'- + ^e^^), (7) 

\ V V w / 

with V = 2- Dh+l3 and w = {1 + X)^{1 + 2X){1 + DhX). 
The value of pc is determined by the definition of the 
front at z = 0, i.e. ^1\"(0) 1. 

Higher order corrections to the analytical profiles 
given above can be computed systematically by rewriting 
Eq. (|4]) in the form, 

pum = p^riz) + - / (z - Zl), (8) 

Ps Jo 

where Gp™(z) = ze-^d{z) (with 6{x) denoting the Hcav- 
isidc function) is the Green's function for the linear part 
of Eq. (j4|). Iteration of Eq. ([8|) yields pu{z) as a power se- 
ries in e. The result, together with pi{z) given by Eq. ([3l). 
can then be fed into Eq. ([5]) to give h{z). 

We have carried out the above procedure to the first 
order in e. Figure 2 shows typical h- and p-profiles as 
obtained from our zeroth order (thin solid lines) and first 
order (dashed lines) analytical solution for ps = 4. As 
anticipated earlier, the p-profiles (black) for the motile 
population have the shape of a bulge with a depletion 
zone right ahead of the front at z = 0. In the zeroth 
order approximation, the bulge is located at z = 1 with 
a peak value pj^ = 2pc/e ~ 0. 736 pc- For the values of 
Dh and $ shown, the AHL profiles (red) also develop a 
dip in the contact zone. Nonetheless, the traveling wave 
solutions are self-consistent as h never cross the threshold 
(dotted horizontal line) on the motile side. 
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FIG. 2: (Color online.) Profiles of scaled cell density p{z) 
(black) and AHL concentration h{z) (red) around the edge 
of the advancing aggregate at S = 0. Shown here are the 
analytical solution to the zeroth order (thin solid lines) and 
first order (dashed lines) in e — Pc/ps, and the numerically 
exact solution to the steady traveling-wave equations (thick 
solid lines). Here Dh = Du/Dp = 1, /3 = /3/7 = 4, and 
ps = aps/Pho = 4. 

To test our analytical solution we have calculated the 
steady traveling profiles by integrating Eqs. ([I])-© nu- 
merically in the moving frame for the above boundary 
conditions (thick solid lines). As is evident from Fig. 2, 
the zeroth order solution already captures the key fea- 
tures of the solution while the first order solution shows 
quantitatively excellent agreement even at pg = 4. 

Given this good agreement, we can now use the analyt- 
ical expressions to find the stability limit of the traveling 
wave solution, i.e., parameter values for which the peak 
height hm of ft.ii(z) reaches the threshold value = 1. 
Let us first consider Dh = D^/Dp ~ 1 as in the ex- 
periments. The Green's function ^ decays at a rate 
of order one in scale d unit s when the scaled AHL diffu- 
sion length Lfi = \J Dh/ P ^ ^ ~ Lp, but much faster 

when ^ 1. In the latter case, the AHL profile follows 
closely the cell density profile, reaching its peak value at 
the tip of the bulge. A straightforward exercise based on 
Eq. ([7]) of the linear case shows hm = Pm = '^Pc/e while 
h{{)) = pc/2 = 1. Hence hm = 4/e ~ 1.47 > ho- In this 
case the traveling wave solution is not self-consistent. An 
increase of Lh allows immotile cells to contribute more to 
the AHL level in the motile region. Consequently hu{z) 
flattens while pc decreases at the same time. Eventually 
hm drops to a value below Kh to restore self-consistency 
of the traveling wave solution. 

The actual stability limit can be obtained by numer- 
ically solving the equations /iii(zm) = hm = 1 and 
dzhu{zm) = 0, where Zm is the peak position of the 
AHL profile. Using the respective analytical profiles, 
we obtain the zeroth order (3 = (f)^™{Dh) (thin solid 
line) and first order /3 = (j)^^^Dh) (dashed line) phase 
boundaries as shown in Fig. [3l In the latter case, the 
first order AHL profile allows us to compute the shift 
6/3 = —etp{Dh)(l>^^'^{Dh) in /3 that satisfies these equations 
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FIG. 3: Phase diagram for stripe formation. The thin solid 
and dashed curves are, respectively, the phase boundaries as 
calculated from our analytical solution in zeroth and first or- 
der in e = pel ps- The thick solid line is obtained from the 
numerical solution of Eqs. (HJ-lIJ]) in the moving frame. The 
red line represents the boundary determined from the onset of 
stripe patterns based on simulation of the full kinetic model 
([T])-©. Inset: thin and thick solid lines give Lh/Lp on the 
zeroth order and exact phase boundaries, respectively. The 
dotted line shows the function ^{Dh/ Dp) from the first order 
correction in e to the boundary position. 



to order e at a given Dh- The modified boundary is then 
obtained from <i)^^\]Dh) = <?!)""(i)/i) exp[-e?/'(i)/»)]. The 
function ip(Dh) is given by the dotted line in the inset 
of Fig. 3. As a comparison, we have also computed the 
phase boundary /3 — (l){Dh) where hm = 1 using the nu- 
merically exact traveling wave solution (thick solid line 
in Fig. 3). The agreement with the first order phase 
boundary 4>^^\Dh) is very good. 

As a confirmation that our ansatz indeed captures the 
dynamic instability behind the stripe formation process, 
we also show in Fig. 3 (red dots) the actual onset of 
stripes observed from a numerical simulation of Eqs. ([1])- 
([2]). Due to the time it takes for transient stripes to 
dissipate close to the transition with the setup of Fig. 
[IJb), the simulation tends to underestimate the no-stripe 
region. Thus the true phase boundary in the long-time 
limit is expected to be somewhat above the red line. 

This study has led to the following picture of the stripe 
formation process: the growth and lateral expansion of 
the colony into an initially cell free region is described 
by a traveling wave solution. In the steadily propagating 
case, the density-coupled cell motility control breaks the 
colony into an immotile region behind a moving bound- 
ary and a density bulge of motile cells ahead of it. Max- 
imum cell density in the motile region is reached at a 
distance Lp = -v/ DqH from the boundary. In the exper- 
iments of Ref. [l5| , the density coupling is implemented 
via a small molecule AHL which provides information on 



cell density within a distance Lh = \/~DhT^- We have 
shown that the steadily propagating wave is stable when 
L/j is greater than or comparable to Lp. In the oppo- 
site case Lh Lp^ instability develops as the maximum 
AHL concentration in the motile region would exceed the 
threshold /iq for motility suppression. Instead, the colony 
expands with periodic nucleation of new immotile regions 
within the motile bulge ahead of the previously formed 
high-density strip. Cell density behind the moving front 
continue to grow until nutrient exhaustion, where the 
density modulation becomes frozen. From the inset of 
Fig. 3 wc see that the ratio Lh/Lp generally lies around 
0.5 on the phase boundary between the two regimes. 

In our system the propagating front thus drives sequen- 
tial stripe formation in an open geometry. This is very 
different from the classical Swift-Hohenberg [l^ mecha- 
nism where finite-wavelength symmetry breaking insta- 
bility develops in the bulk. The highly nonlinear and 
localized process in the nucleation of new stripes also 
makes our mechanism different from that of pattern for- 
mation driven by fronts propagating into a bistable sys- 
tem where modulations arise during the linear instability 
development at the front [2^. In this respect, there are 
some similarities between our system and the nonperi- 
odic Liesegang patterns since in both cases new "phase" 



precipitates when certain critical density is reached [21 



On the other hand, in the chemical systems that exhibit 
Liesegang patterns, reactant density increases via trans- 
port instead of growth. 
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Supporting Material: Stripe formation in bacterial systems with 

density-dependent motility 

As mentioned our model can be derived by coarse-graining of the microscopic dynamics. Generally, this procedure 
yields a mixture of two terms dx{iJ,{h)dxp) + Odx{pdxfJ,{h)), where 6 depends on the underlying microscopic dynamics, 
e.g. 6* = 1 for Ito-dynamics and 9 = 1/2 for Stratonovich dynamics With this general coupling our model becomes 



dh 
'dt 
dp 
dt 



+ ap - (3h, 



— i^,ih)p)-{i-e)-ip^ 



7op 1 

Ps 



(SI) 
(S2) 



In the main text, we study the case 9 = 1. In the following, we demonstrate that for 9 < 1 the moving front still 
acts as an absorbing boundary and that our main conclusions remain valid in this case. 
Boundary conditions at the moving front z = 0. 

The field h{x) is continuous at ft, = Kh and even differentiable. This can be seen by integrating Eq. (|S1|) from 
z = 0~ to z = O'^ (as in the main text z = x — ct) and by using dh/dt = —cdh/dz. 



cw = -c[h{0+) - h{0-)] = Dh dh/dz 



z=0+ 

2=0- ' 



implying dh{z)/dz\^+ = dh{z)/dz\f^- as w — > 0. 
In contrast, the density p{z) is discontinuous at z 



with 



(S3) 
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where we have used that p{h{z = 0+)) = p{h{z = 0^)). Eq. ([S4| thus implies /9(0+) = as w — > 0. Thus, 
independent of the value of 9 the moving front acts as absorbing boundary. From this analysis it is also evident 
that the boundary condition p(0"^) = does not depend on our assumptions on the specific functional form of the 
interpolating function Dp{h) for Kh — w < h < Kh- 

Finally, the slope of p(0+) can be determined by integrating Eq. (|Sip from z = 0^ to z = 0+ 



dp/dz 



=0+ 



cpc/D^ 



(S5) 



where (as in the main text) pc = p(0~). 

Results for 9 < 1. We have calculated the cellular density profile and the phase diagram for different values 9 < 1. 
As can be seen from Fig. [ST] the additional term in Eq. (|S2p only leads to very small (hardly visible) modifications 
of the cellular density and AHL concentration profile. Consequently, the phase boundary is also not affected by the 
value of 9. 
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FIG. SI: Profiles of scaled cellular density p{z) (black) and AHL concentration h[z) (red) close to the moving front at S = 
as calculated numerically from Eqs. (|S1[) and (|S2|l . The solid lines are for 6 = 1/2 the dashed lines for 6 = 1. Data are for 
Dh = Dh/Dp = 1.3, j3 = Ph = 5, and p» = aps/Pho = 4. 



